clear all;

% Keyword 1

n=6;
m=0.31;
v=.26^2;
mu = log((m^2)/sqrt(v+m^2));
sigma = sqrt(log(v/(m^2)+1));
fun = @(x)(-max(0.00001,x)*(1-normcdf(log(max(0.00001,x)),mu,sigma)));
x0 = 1;
rstar = fminsearch(fun,x0)

nobs = 10000000; 
vals = lognrnd(mu,sigma,nobs,n);
vals = sort(vals,2);

alpha = .7;
refctr = [1];
for i = 1:n-1
    refctr = [refctr;alpha^i];
end;
    
aux = repmat(sort(refctr'),nobs,1);

refctr2 = refctr.*(1:n)';
refctr2 = (1-alpha)*refctr2(1:n-1);
aux2 = repmat(fliplr(refctr2'),nobs,1);

r=0;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e0 = mean(sum(revsdirect,2)+aux3)

r=.1;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e10 = mean(sum(revsdirect,2)+aux3)

r=(.1+rstar)/2;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eMP = mean(sum(revsdirect,2)+aux3)

r=rstar;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eOPT = mean(sum(revsdirect,2)+aux3)

% Keyword 2

n=6;
m=0.78;
v=.11^2;
mu = log((m^2)/sqrt(v+m^2));
sigma = sqrt(log(v/(m^2)+1));
fun = @(x)(-max(0.00001,x)*(1-normcdf(log(max(0.00001,x)),mu,sigma)));
x0 = 1;
rstar = fminsearch(fun,x0)

nobs = 10000000; 
vals = lognrnd(mu,sigma,nobs,n);
vals = sort(vals,2);

alpha = .7;
refctr = [1];
for i = 1:n-1
    refctr = [refctr;alpha^i];
end;
    
aux = repmat(sort(refctr'),nobs,1);

refctr2 = refctr.*(1:n)';
refctr2 = (1-alpha)*refctr2(1:n-1);
aux2 = repmat(fliplr(refctr2'),nobs,1);

r=0;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e0 = mean(sum(revsdirect,2)+aux3)

r=.1;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e10 = mean(sum(revsdirect,2)+aux3)

r=(.1+rstar)/2;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eMP = mean(sum(revsdirect,2)+aux3)

r=rstar;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eOPT = mean(sum(revsdirect,2)+aux3)


% Keyword 3

n=5;
m=0.22;
v=.13^2;
mu = log((m^2)/sqrt(v+m^2));
sigma = sqrt(log(v/(m^2)+1));
fun = @(x)(-max(0.00001,x)*(1-normcdf(log(max(0.00001,x)),mu,sigma)));
x0 = 1;
rstar = fminsearch(fun,x0)

nobs = 10000000; 
vals = lognrnd(mu,sigma,nobs,n);
vals = sort(vals,2);

alpha = .7;
refctr = [1];
for i = 1:n-1
    refctr = [refctr;alpha^i];
end;
    
aux = repmat(sort(refctr'),nobs,1);

refctr2 = refctr.*(1:n)';
refctr2 = (1-alpha)*refctr2(1:n-1);
aux2 = repmat(fliplr(refctr2'),nobs,1);

r=0;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e0 = mean(sum(revsdirect,2)+aux3)

r=.1;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e10 = mean(sum(revsdirect,2)+aux3)

r=(.1+rstar)/2;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eMP = mean(sum(revsdirect,2)+aux3)

r=rstar;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eOPT = mean(sum(revsdirect,2)+aux3)

% Keyword 4

n=9;
m=2.18;
v=.99^2;
mu = log((m^2)/sqrt(v+m^2));
sigma = sqrt(log(v/(m^2)+1));
fun = @(x)(-max(0.00001,x)*(1-normcdf(log(max(0.00001,x)),mu,sigma)));
x0 = 1;
rstar = fminsearch(fun,x0)

nobs = 10000000; 
vals = lognrnd(mu,sigma,nobs,n);
vals = sort(vals,2);

alpha = .7;
refctr = [1];
for i = 1:n-1
    refctr = [refctr;alpha^i];
end;
    
aux = repmat(sort(refctr'),nobs,1);

refctr2 = refctr.*(1:n)';
refctr2 = (1-alpha)*refctr2(1:n-1);
aux2 = repmat(fliplr(refctr2'),nobs,1);

r=0;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e0 = mean(sum(revsdirect,2)+aux3)

r=.1;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e10 = mean(sum(revsdirect,2)+aux3)

r=(.1+rstar)/2;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eMP = mean(sum(revsdirect,2)+aux3)

r=rstar;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eOPT = mean(sum(revsdirect,2)+aux3)